/*
 Copyright (c) 2011, The Cinder Project, All rights reserved.
 This code is intended for use with the Cinder C++ library: http://libcinder.org

 Portions Copyright (c) 2004, Laminar Research.

 Redistribution and use in source and binary forms, with or without modification, are permitted provided that
 the following conditions are met:

    * Redistributions of source code must retain the above copyright notice, this list of conditions and
	the following disclaimer.
    * Redistributions in binary form must reproduce the above copyright notice, this list of conditions and
	the following disclaimer in the documentation and/or other materials provided with the distribution.

 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED
 WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
 PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR
 ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
 TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
 HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
 NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 POSSIBILITY OF SUCH DAMAGE.

 Portions Copyright Geometric Tools LLC, Redmond WA 98052
 Boost Software License - Version 1.0 - August 17th, 2003

 Permission is hereby granted, free of charge, to any person or organization
 obtaining a copy of the software and accompanying documentation covered by
 this license (the "Software") to use, reproduce, display, distribute,
 execute, and transmit the Software, and to prepare derivative works of the
 Software, and to permit third-parties to whom the Software is furnished to
 do so, all subject to the following:

 The copyright notices in the Software and this entire statement, including
 the above license grant, this restriction and the following disclaimer,
 must be included in all copies of the Software, in whole or in part, and
 all derivative works of the Software, unless such copies or derivative
 works are solely in the form of machine-executable object code generated by
 a source language processor.

 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 FITNESS FOR A PARTICULAR PURPOSE, TITLE AND NON-INFRINGEMENT. IN NO EVENT
 SHALL THE COPYRIGHT HOLDERS OR ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE
 FOR ANY DAMAGES OR OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE,
 ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
 DEALINGS IN THE SOFTWARE.
*/

#include "cinder/CinderMath.h"
#include <algorithm>
#include <vector>

using namespace glm;

namespace cinder {


/////////////////////////////////////////////////////////////////////////////////////////////////
// solveCubic
template<typename T>
int solveCubic( T a, T b, T c, T d, T result[3] )
{
	if( a == 0 )
		return solveQuadratic( b, c, d, result );

	T a2 = a * a;
	T b2 = b * b;
	T f = ((3 * c / a) - (b2 / a2)) / 3;
	T g = ((2 * b2 * b) / (a2 * a) - (9 * b * c) / a2 + (27 * d) / a) / 27;
	T h = g * g / 4 + f * f * f / 27;

	if( f == 0 && g == 0 && h == 0 ) {
		result[0] = -math<T>::cbrt( d / a );
		return 1;
	}
	else if( h > 0 ) {
		// 1 root
		T hr = math<T>::sqrt( h );
		T r = -( g / 2 ) + hr;
		T s = math<T>::cbrt( r );
		T t = -(g / 2) - hr;
		T u = math<T>::cbrt( t );

		result[0] = (s + u) - (b / (3 * a));
		return 1;
	}
	else { // 3 roots
		T i = math<T>::sqrt( (g * g / 4) - h );
		T j = math<T>::cbrt( i );
		T k = math<T>::acos( -(g / (2 * i)) );
		T l = -j;
		T m = math<T>::cos( k / 3 );
		T n = math<T>::sqrt(3) * math<T>::sin( k / 3 );
		T p = -b / (3 * a);
		result[0] = 2 * j * m - (b / (3 * a));
		result[1] = l * (m + n) + p;
		result[2] = l * (m - n) + p;
		return 3;
	}
}
template CI_API int solveCubic( float a, float b, float c, float d, float result[3] );
template CI_API int solveCubic( double a, double b, double c, double d, double result[3] );

namespace {
float PointOnEllipseBisector( int numComponents, const vec2 &extents, const vec2 &y, vec2& x )
{
	vec2 z;
	float sumZSqr = 0;
	int i;
	for( i = 0; i < numComponents; ++i ) {
		z[i] = y[i] / extents[i];
		sumZSqr += z[i] * z[i];
	}

	if( sumZSqr == 1 ) {
		// The point is on the hyperellipsoid.
		for (i = 0; i < numComponents; ++i)
			x[i] = y[i];

		return 0;
	}

	float emin = extents[numComponents - 1];
	vec2 pSqr, numerator;
	for( i = 0; i < numComponents; ++i ) {
		float p = extents[i] / emin;
		pSqr[i] = p * p;
		numerator[i] = pSqr[i] * z[i];
	}

	// The maximum number of bisections required for Real before the interval
	// endpoints are equal (as floating-point numbers).
	int const jmax = std::numeric_limits<float>::digits - std::numeric_limits<float>::min_exponent;

	float s = 0, smin = z[numComponents - 1] - 1, smax;
	if( sumZSqr < 1 )
		// The point is strictly inside the hyperellipsoid.
		smax = 0;
	else
		// The point is strictly outside the hyperellipsoid.
		//smax = LengthRobust(numerator) - (Real)1;
		smax = length( numerator ) - 1;

	for( int j = 0; j < jmax; ++j ) {
		s = (smin + smax) * 0.5f;
		if (s == smin || s == smax)
			break;

		float g = -1;
		for( i = 0; i < numComponents; ++i ) {
			float ratio = numerator[i] / (s + pSqr[i]);
			g += ratio * ratio;
		}

		if( g > 0 )
			smin = s;
		else if( g < 0 )
			smax = s;
		else
			break;
	}

	float sqrDistance = 0;
	for( i = 0; i < numComponents; ++i ) {
		x[i] = pSqr[i] * y[i] / (s + pSqr[i]);
		float diff = x[i] - y[i];
		sqrDistance += diff*diff;
	}
	return sqrDistance;
}

float PointOnEllipseSqrDistanceSpecial( const vec2 &extents, const vec2 &y, vec2 &x )
{
    float sqrDistance = 0;

    vec2 ePos, yPos, xPos;
    int numPos = 0;
    for( int i = 0; i < 2; ++i ) {
        if( y[i] > 0 ) {
            ePos[numPos] = extents[i];
            yPos[numPos] = y[i];
            ++numPos;
        }
        else
            x[i] = 0;
    }

    if( y[2 - 1] > 0 )
        sqrDistance = PointOnEllipseBisector( numPos, ePos, yPos, xPos );
    else {  // y[N-1] = 0
        float numer[1], denom[1];
        float eNm1Sqr = extents[2 - 1] * extents[2 - 1];
        for( int i = 0; i < numPos; ++i)
        {
            numer[i] = ePos[i] * yPos[i];
            denom[i] = ePos[i] * ePos[i] - eNm1Sqr;
        }

        bool inSubHyperbox = true;
        for( int i = 0; i < numPos; ++i) {
            if( numer[i] >= denom[i]) {
                inSubHyperbox = false;
                break;
            }
        }

        bool inSubHyperellipsoid = false;
        if( inSubHyperbox ) {
            // yPos[] is inside the axis-aligned bounding box of the
            // subhyperellipsoid.  This intermediate test is designed to guard
            // against the division by zero when ePos[i] == e[N-1] for some i.
            float xde[1];
            float discr = 1;
            for( int i = 0; i < numPos; ++i)
            {
                xde[i] = numer[i] / denom[i];
                discr -= xde[i] * xde[i];
            }
            if( discr > 0 ) {
                // yPos[] is inside the subhyperellipsoid.  The closest
                // hyperellipsoid point has x[N-1] > 0.
                sqrDistance = 0;
                for( int i = 0; i < numPos; ++i)
                {
                    xPos[i] = ePos[i] * xde[i];
                    float diff = xPos[i] - yPos[i];
                    sqrDistance += diff*diff;
                }
                x[2 - 1] = extents[2 - 1] * sqrt(discr);
                sqrDistance += x[2 - 1] * x[2 - 1];
                inSubHyperellipsoid = true;
            }
        }

        if( ! inSubHyperellipsoid ) {
            // yPos[] is outside the subhyperellipsoid.  The closest
            // hyperellipsoid point has x[N-1] == 0 and is on the
            // domain-boundary hyperellipsoid.
            x[2 - 1] = 0;
            sqrDistance = PointOnEllipseBisector( numPos, ePos, yPos, xPos );
        }
    }

    // Fill in those x[] values that were not zeroed out initially.
    for( int i = 0, numPos = 0; i < 2; ++i ) {
        if( y[i] > 0 ) {
            x[i] = xPos[numPos];
            ++numPos;
        }
    }

    return sqrDistance;
}

float PointOnEllipseSqrDistance( const vec2 &extents, const vec2 &y, vec2 &x )
{
    // Determine negations for y to the first octant.
    bool negate[2];
    for( int i = 0; i < 2; ++i )
        negate[i] = y[i] < 0;

    // Determine the axis order for decreasing extents.
    std::pair<float, int> permute[2];
    for( int i = 0; i < 2; ++i ) {
        permute[i].first = -extents[i];
        permute[i].second = i;
    }
    std::sort( &permute[0], &permute[2] );

    int invPermute[2];
    for( int i = 0; i < 2; ++i )
        invPermute[permute[i].second] = i;

    vec2 locE, locY;
	int j;
    for( int i = 0; i < 2; ++i ) {
        j = permute[i].second;
        locE[i] = extents[j];
        locY[i] = std::abs(y[j]);
    }

    vec2 locX;
    float sqrDistance = PointOnEllipseSqrDistanceSpecial( locE, locY, locX );

    // Restore the axis order and reflections.
    for( int i = 0; i < 2; ++i ) {
        j = invPermute[i];
        if( negate[i] )
            locX[j] = -locX[j];
        x[i] = locX[j];
    }

    return sqrDistance;
}
} // anonymous namespace for closestPointOnEllipse

vec2 getClosestPointEllipse( const vec2& center, const vec2& axisA, const vec2& axisB, const vec2& testPoint )
{
	// Compute the coordinates of Y in the hyperellipsoid coordinate system.
	float lengthA = length( axisA );
	float lengthB = length( axisB );
	vec2 unitA = axisA / lengthA;
	vec2 unitB = axisB / lengthB;
	vec2 diff = testPoint - center;
	vec2 y( dot( diff, unitA ), dot( diff, unitB ) );

	// Compute the closest hyperellipsoid point in the axis-aligned
	// coordinate system.
	vec2 x;
	vec2 extents( lengthA, lengthB );
	PointOnEllipseSqrDistance( extents, y, x );

	// Convert back to the original coordinate system.
	vec2 result = center;
	result += x[0] * unitA;
	result += x[1] * unitB;

	return result;
}

namespace {

//! Given a point and a Bezier curve, generate a 5th-degree Bezier - format equation whose solution finds the point on the curve nearest the user - defined point.
template<typename T>
void convertToBezierForm( const glm::tvec2<T, glm::defaultp> *controlPoints, const glm::tvec2<T, glm::defaultp> &testPoint, glm::tvec2<T, glm::defaultp> *w )
{
	static const T z[3][4] = { { T(1.0), T(0.6), T(0.3), T(0.1) },{ T(0.4), T(0.6), T(0.6), T(0.4) },{ T(0.1), T(0.3), T(0.6), T(1.0) } };

	// Determine the c's -- these are vectors created by subtracting point P from each of the control points.
	glm::tvec2<T, glm::defaultp> c[4];
	c[0] = controlPoints[0] - testPoint;
	c[1] = controlPoints[1] - testPoint;
	c[2] = controlPoints[2] - testPoint;
	c[3] = controlPoints[3] - testPoint;

	// Determine the d's -- these are vectors created by subtracting each control point from the next.
	glm::tvec2<T, glm::defaultp> d[3];
	d[0] = ( controlPoints[1] - controlPoints[0] ) * T{3};
	d[1] = ( controlPoints[2] - controlPoints[1] ) * T{3};
	d[2] = ( controlPoints[3] - controlPoints[2] ) * T{3};

	// Create the c,d table -- this is a table of dot products of the c's and d's.
	T cdTable[3][4];
	for( int row = 0; row <= 2; row++ ) {
		for( int column = 0; column <= 3; column++ ) {
			cdTable[row][column] = glm::dot( d[row], c[column] );
		}
	}

	// Now, apply the z's to the dot products, on the skew diagonal. Also, set up the x-values, making these "points".
	for( int i = 0; i <= 5; i++ ) {
		w[i].y = 0.0;
		w[i].x = (T)( i ) / 5;
	}

	const int n = 3;
	const int m = 2;
	for( int k = 0; k <= n + m; k++ ) {
		int lb = glm::max( 0, k - m );
		int ub = glm::min( k, n );
		for( int i = lb; i <= ub; i++ ) {
			int j = k - i;
			w[k].y += cdTable[j][i] * z[j][i];
		}
	}
}

//! Returns the sign of the value (-1 if < 0, 0 if == 0, 1 if > 0). Faster than using glm::sign()!
template <typename T> 
int sgn( T val )
{
	return ( T( 0 ) < val ) - ( val < T( 0 ) );
}

//! Count the number of times a Bezier control polygon crosses the x-axis. This number is >= the number of roots.
template<typename T>
int crossingCount( const glm::tvec2<T, glm::defaultp> *controlPoints, int degree )
{
	int n = 0;

	int sign = sgn( controlPoints[0].y );
	for( int i = 1; i <= degree; i++ ) {
		int s = sgn( controlPoints[i].y );
		if( sign != s )
			n++;
		sign = s;
	}

	return n;
}

//! Check if the control polygon of a Bezier curve is flat enough for recursive subdivision to bottom out.
template<typename T>
int controlPolygonFlatEnough( const glm::tvec2<T, glm::defaultp> *controlPoints, int degree )
{
	// Derive the implicit equation for line connecting first and last control points.
	T a = controlPoints[0].y - controlPoints[degree].y;
	T b = controlPoints[degree].x - controlPoints[0].x;
	T c = controlPoints[0].x * controlPoints[degree].y - controlPoints[degree].x * controlPoints[0].y;

	T max_distance_above{0};
	T max_distance_below{0};

	T value;
	for( int i = 1; i < degree; i++ ) {
		value = a * controlPoints[i].x + b * controlPoints[i].y + c;
		max_distance_above = glm::max( max_distance_above, value );
		max_distance_below = glm::min( max_distance_below, value );
	}

	// Compute intercepts of bounding box.
	T error = glm::abs( ( max_distance_above - max_distance_below ) / a );

	static const T kEpsilon = std::numeric_limits<T>::epsilon();
	return ( error < kEpsilon ) ? 1 : 0;
}

//! Compute intersection of chord from first control point to last with x-axis.
template<typename T>
T computeXIntercept( const glm::tvec2<T, glm::defaultp> &p0, const glm::tvec2<T, glm::defaultp> &p3 )
{
	glm::tvec2<T, glm::defaultp> d = p3 - p0;
	return ( d.x * p0.y - d.y * p0.x ) / ( -d.y );
}


// Evaluate a Bezier curve at a particular parameter value. Assumes the curve is cubic (degree == 3).
template<typename T>
glm::tvec2<T, glm::defaultp> bezier( const glm::tvec2<T, glm::defaultp> *controlPoints, T t )
{
	T u = T{ 1 } -t;
	return u * u * ( u * controlPoints[0] + T{ 3 } * t * controlPoints[1] ) + t * t * ( T{ 3 } * u * controlPoints[2] + t * controlPoints[3] );
}

// Evaluate a Bezier curve at a particular parameter value. Fill in control points for resulting sub-curves in \a left and \a right.
template<typename T>
glm::tvec2<T, glm::defaultp> bezier( const glm::tvec2<T, glm::defaultp> *controlPoints, int degree, T t, glm::tvec2<T, glm::defaultp> *left, glm::tvec2<T, glm::defaultp> *right )
{
	assert( nullptr != left );
	assert( nullptr != right );
	assert( degree <= 5 );

	glm::tvec2<T, glm::defaultp> value[6];

	// Copy control points.
	memcpy( value, controlPoints, sizeof( value ) );

	// Initialize the first sub-curve.
	left[0] = controlPoints[0];

	// Triangle computation.
	for( int i = 1; i <= degree; i++ ) {
		for( int j = 0; j <= degree - i; j++ ) {
			value[j] += t * ( value[j + 1] - value[j] );
		}

		// Copy parameters for the first sub-curve.
		left[i] = value[0];
	}

	// Copy parameters for the second sub-curve.
	memcpy( right, value, sizeof( value ) );

	return value[degree];
}

//! Given a 5th-degree equation in Bernstein-Bezier form, find all of the roots in the interval[0, 1]. Return the number of roots found.
template<typename T>
int findRoots( const glm::tvec2<T, glm::defaultp> *controlPoints, int degree, T *t, int depth )
{
	switch( crossingCount( controlPoints, degree ) ) {
		case 0: { // No solutions here.
			return 0;
		}
		case 1: { // Unique solution
			// Stop recursion when the tree is deep enough. If deep enough, return 1 solution at midpoint.
			const int kMaxDepth = 64;
			if( depth >= kMaxDepth ) {
				t[0] = ( controlPoints[0].x + controlPoints[5].x ) / T{2};
				return 1;
			}
			if( controlPolygonFlatEnough<T>( controlPoints, degree ) ) {
				t[0] = computeXIntercept<T>( controlPoints[0], controlPoints[degree] );
				return 1;
			}
			break;
		}
	}

	// Otherwise, solve recursively after subdividing control polygon.
	glm::tvec2<T, glm::defaultp> left[6], right[6];
	bezier<T>( controlPoints, degree, 0.5f, left, right );
	
	T left_t[6], right_t[6];
	int left_count = findRoots<T>( left, degree, left_t, depth + 1 );
	int right_count = findRoots<T>( right, degree, right_t, depth + 1 );

	// Gather solutions together.
	for( int i = 0; i < left_count; i++ ) {
		t[i] = left_t[i];
	}
	for( int i = 0; i < right_count; i++ ) {
		t[i + left_count] = right_t[i];
	}

	// Send back total number of solutions.
	return ( left_count + right_count );
}

} // anonymous namespace for closestPointCubic

template<typename T>
glm::tvec2<T, glm::defaultp> getClosestPointLinear( const glm::tvec2<T, glm::defaultp> *controlPoints, const glm::tvec2<T, glm::defaultp> & testPoint )
{
	T d = glm::distance2( controlPoints[1], controlPoints[0] );
	if( d > T{0} ) {
		T t = glm::clamp( glm::dot( testPoint - controlPoints[0], controlPoints[1] - controlPoints[0] ) / d, T{0}, T{1} );
		return controlPoints[0] + t * ( controlPoints[1] - controlPoints[0] );
	}
	else {
		return controlPoints[0];
	}
}
template CI_API glm::tvec2<float, glm::defaultp> getClosestPointLinear( const glm::tvec2<float, glm::defaultp> *controlPoints, const glm::tvec2<float, glm::defaultp> & testPoint );
template CI_API glm::tvec2<double, glm::defaultp> getClosestPointLinear( const glm::tvec2<double, glm::defaultp> *controlPoints, const glm::tvec2<double, glm::defaultp> & testPoint );

template<typename T>
glm::tvec2<T, glm::defaultp> getClosestPointQuadratic( const glm::tvec2<T, glm::defaultp> *controlPoints, const glm::tvec2<T, glm::defaultp> & testPoint )
{
	glm::tvec2<T, glm::defaultp> closest;

	glm::tvec2<T, glm::defaultp> M = controlPoints[0] - testPoint;
	glm::tvec2<T, glm::defaultp> A = controlPoints[1] - controlPoints[0];
	glm::tvec2<T, glm::defaultp> B = controlPoints[2] - controlPoints[1] - A;
	T a = glm::dot( B, B );
	T b = 3 * glm::dot( A, B );
	T c = 2 * glm::dot( A, A ) + glm::dot( M, B );
	T d = glm::dot( M, A );
	T t[3];
	int solutions = solveCubic<T>( a, b, c, d, t );

	T distance2 = std::numeric_limits<T>::max();
	for( int i = 0; i < solutions; ++i ) {
		T u = glm::clamp( t[i], T{0}, T{1} );
		T v = T{1} - u;
		glm::tvec2<T, glm::defaultp> p = controlPoints[0] * ( v*v ) + controlPoints[1] * ( T{2} * u*v ) + controlPoints[2] * ( u*u );

		T d = glm::distance2( testPoint, p );
		if( d < distance2 ) {
			closest = p;
			distance2 = d;
		}
	}

	return closest;
}
template CI_API glm::tvec2<float, glm::defaultp> getClosestPointQuadratic( const glm::tvec2<float, glm::defaultp> *controlPoints, const glm::tvec2<float, glm::defaultp> & testPoint );
template CI_API glm::tvec2<double, glm::defaultp> getClosestPointQuadratic( const glm::tvec2<double, glm::defaultp> *controlPoints, const glm::tvec2<double, glm::defaultp> & testPoint );

template<typename T>
glm::tvec2<T, glm::defaultp> getClosestPointCubic( const glm::tvec2<T, glm::defaultp> *controlPoints, const glm::tvec2<T, glm::defaultp> & testPoint )
{
	// Convert problem to 5th-degree Bezier form.
	glm::tvec2<T, glm::defaultp> w[6];
	convertToBezierForm<T>( controlPoints, testPoint, w );

	// Find all possible roots of 5th-degree equation.
	T t_candidate[5];
	int n_solutions = findRoots<T>( w, 5, t_candidate, 0 );

	// Compare distances of P to all candidates, and to t=0, and t=1.
	T t{ 0 };
	{
		// Check distance to beginning of curve, where t = 0.
		T dist = glm::distance2( testPoint, controlPoints[0] );

		// Find distances for candidate points.
		glm::tvec2<T, glm::defaultp> p;
		for( int i = 0; i < n_solutions; i++ ) {
			p = bezier<T>( controlPoints, t_candidate[i] );
			T new_dist = glm::distance2( testPoint, p );
			if( new_dist < dist ) {
				dist = new_dist;
				t = t_candidate[i];
			}
		}

		// Finally, look at distance to end point, where t = 1.0.
		T new_dist = glm::distance2( testPoint, controlPoints[3] );
		if( new_dist < dist ) {
			dist = new_dist;
			t = T{ 1 };
		}
	}

	// Return the point on the curve at parameter value t.
	return ( bezier<T>( controlPoints, t ) );
}
template CI_API glm::tvec2<float, glm::defaultp> getClosestPointCubic<float>( const glm::tvec2<float, glm::defaultp> *controlPoints, const glm::tvec2<float, glm::defaultp> & testPoint );
template CI_API glm::tvec2<double, glm::defaultp> getClosestPointCubic<double>( const glm::tvec2<double, glm::defaultp> *controlPoints, const glm::tvec2<double, glm::defaultp> & testPoint );

union float32_t
{
	float f;
	uint u;
	struct {
		uint Mantissa : 23;
		uint Exponent : 8;
		uint Sign : 1;
	};
};

// Algorithm due to Fabian "ryg" Giesen.
static half_float float_to_half( float32_t f )
{
    float32_t f32infty = { 255 << 23 };
    float32_t f16infty = { 31 << 23 };
    float32_t magic = { 15 << 23 };
    uint sign_mask = 0x80000000u;
    uint round_mask = ~0xfffu; 
    half_float o = { 0 };
 
    uint sign = f.u & sign_mask;
    f.u ^= sign;
 
    // NOTE all the integer compares in this function can be safely
    // compiled into signed compares since all operands are below
    // 0x80000000. Important if you want fast straight SSE2 code
    // (since there's no unsigned PCMPGTD).
 
    if (f.u >= f32infty.u) // Inf or NaN (all exponent bits set)
        o.u = (f.u > f32infty.u) ? 0x7e00 : 0x7c00; // NaN->qNaN and Inf->Inf
    else // (De)normalized number or zero
    {
        f.u &= round_mask;
        f.f *= magic.f;
        f.u -= round_mask;
        if (f.u > f16infty.u) f.u = f16infty.u; // Clamp to signed infinity if overflowed
 
        o.u = f.u >> 13; // Take the bits!
    }
 
    o.u |= sign >> 16;
    return o;
}

cinder::half_float floatToHalf( float f )
{
	return float_to_half( float32_t( { f } ) );
}

// Algorithm due to Fabian "ryg" Giesen.
float halfToFloat( cinder::half_float h )
{
	static const float32_t magic = { 113 << 23 };
	static const uint shifted_exp = 0x7c00 << 13; // exponent mask after shift
	float32_t o;

	o.u = (h.u & 0x7fff) << 13;     // exponent/mantissa bits
	uint exp = shifted_exp & o.u;   // just the exponent
	o.u += (127 - 15) << 23;        // exponent adjust

	// handle exponent special cases
	if (exp == shifted_exp) // Inf/NaN?
		o.u += (128 - 16) << 23;    // extra exp adjust
	else if (exp == 0) { // Zero/Denormal?
		o.u += 1 << 23;             // extra exp adjust
		o.f -= magic.f;             // renormalize
	}

	o.u |= (h.u & 0x8000) << 16;    // sign bit
	return o.f;
}

} // namespace cinder